home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / quadpack / dqagi.f next >
Text File  |  1996-07-19  |  9KB  |  192 lines
  1.       SUBROUTINE DQAGI(F,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,
  2.      *   IER,LIMIT,LENW,LAST,IWORK,WORK)
  3. C***BEGIN PROLOGUE  DQAGI
  4. C***DATE WRITTEN   800101   (YYMMDD)
  5. C***REVISION DATE  830518   (YYMMDD)
  6. C***CATEGORY NO.  H2A3A1,H2A4A1
  7. C***KEYWORDS  AUTOMATIC INTEGRATOR, INFINITE INTERVALS,
  8. C             GENERAL-PURPOSE, TRANSFORMATION, EXTRAPOLATION,
  9. C             GLOBALLY ADAPTIVE
  10. C***AUTHOR  PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
  11. C           DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. -K.U.LEUVEN
  12. C***PURPOSE  THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  13. C            INTEGRAL   I = INTEGRAL OF F OVER (BOUND,+INFINITY)
  14. C            OR I = INTEGRAL OF F OVER (-INFINITY,BOUND)
  15. C            OR I = INTEGRAL OF F OVER (-INFINITY,+INFINITY)
  16. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  17. C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  18. C***DESCRIPTION
  19. C
  20. C        INTEGRATION OVER INFINITE INTERVALS
  21. C        STANDARD FORTRAN SUBROUTINE
  22. C
  23. C        PARAMETERS
  24. C         ON ENTRY
  25. C            F      - DOUBLE PRECISION
  26. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  27. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  28. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  29. C
  30. C            BOUND  - DOUBLE PRECISION
  31. C                     FINITE BOUND OF INTEGRATION RANGE
  32. C                     (HAS NO MEANING IF INTERVAL IS DOUBLY-INFINITE)
  33. C
  34. C            INF    - INTEGER
  35. C                     INDICATING THE KIND OF INTEGRATION RANGE INVOLVED
  36. C                     INF = 1 CORRESPONDS TO  (BOUND,+INFINITY),
  37. C                     INF = -1            TO  (-INFINITY,BOUND),
  38. C                     INF = 2             TO (-INFINITY,+INFINITY).
  39. C
  40. C            EPSABS - DOUBLE PRECISION
  41. C                     ABSOLUTE ACCURACY REQUESTED
  42. C            EPSREL - DOUBLE PRECISION
  43. C                     RELATIVE ACCURACY REQUESTED
  44. C                     IF  EPSABS.LE.0
  45. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  46. C                     THE ROUTINE WILL END WITH IER = 6.
  47. C
  48. C
  49. C         ON RETURN
  50. C            RESULT - DOUBLE PRECISION
  51. C                     APPROXIMATION TO THE INTEGRAL
  52. C
  53. C            ABSERR - DOUBLE PRECISION
  54. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  55. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  56. C
  57. C            NEVAL  - INTEGER
  58. C                     NUMBER OF INTEGRAND EVALUATIONS
  59. C
  60. C            IER    - INTEGER
  61. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  62. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  63. C                             ACCURACY HAS BEEN ACHIEVED.
  64. C                   - IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE. THE
  65. C                             ESTIMATES FOR RESULT AND ERROR ARE LESS
  66. C                             RELIABLE. IT IS ASSUMED THAT THE REQUESTED
  67. C                             ACCURACY HAS NOT BEEN ACHIEVED.
  68. C            ERROR MESSAGES
  69. C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  70. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
  71. C                             SUBDIVISIONS BY INCREASING THE VALUE OF
  72. C                             LIMIT (AND TAKING THE ACCORDING DIMENSION
  73. C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
  74. C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
  75. C                             TO ANALYZE THE INTEGRAND IN ORDER TO
  76. C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
  77. C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
  78. C                             DETERMINED (E.G. SINGULARITY,
  79. C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
  80. C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
  81. C                             INTERVAL AT THIS POINT AND CALLING THE
  82. C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
  83. C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  84. C                             SHOULD BE USED, WHICH IS DESIGNED FOR
  85. C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
  86. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
  87. C                             DETECTED, WHICH PREVENTS THE REQUESTED
  88. C                             TOLERANCE FROM BEING ACHIEVED.
  89. C                             THE ERROR MAY BE UNDER-ESTIMATED.
  90. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
  91. C                             AT SOME POINTS OF THE INTEGRATION
  92. C                             INTERVAL.
  93. C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
  94. C                             ROUNDOFF ERROR IS DETECTED IN THE
  95. C                             EXTRAPOLATION TABLE.
  96. C                             IT IS ASSUMED THAT THE REQUESTED TOLERANCE
  97. C                             CANNOT BE ACHIEVED, AND THAT THE RETURNED
  98. C                             RESULT IS THE BEST WHICH CAN BE OBTAINED.
  99. C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
  100. C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
  101. C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
  102. C                             OF IER.
  103. C                         = 6 THE INPUT IS INVALID, BECAUSE
  104. C                             (EPSABS.LE.0 AND
  105. C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
  106. C                              OR LIMIT.LT.1 OR LENIW.LT.LIMIT*4.
  107. C                             RESULT, ABSERR, NEVAL, LAST ARE SET TO
  108. C                             ZERO. EXEPT WHEN LIMIT OR LENIW IS
  109. C                             INVALID, IWORK(1), WORK(LIMIT*2+1) AND
  110. C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
  111. C                             IS SET TO A AND WORK(LIMIT+1) TO B.
  112. C
  113. C         DIMENSIONING PARAMETERS
  114. C            LIMIT - INTEGER
  115. C                    DIMENSIONING PARAMETER FOR IWORK
  116. C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
  117. C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
  118. C                    (A,B), LIMIT.GE.1.
  119. C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
  120. C
  121. C            LENW  - INTEGER
  122. C                    DIMENSIONING PARAMETER FOR WORK
  123. C                    LENW MUST BE AT LEAST LIMIT*4.
  124. C                    IF LENW.LT.LIMIT*4, THE ROUTINE WILL END
  125. C                    WITH IER = 6.
  126. C
  127. C            LAST  - INTEGER
  128. C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
  129. C                    PRODUCED IN THE SUBDIVISION PROCESS, WHICH
  130. C                    DETERMINES THE NUMBER OF SIGNIFICANT ELEMENTS
  131. C                    ACTUALLY IN THE WORK ARRAYS.
  132. C
  133. C         WORK ARRAYS
  134. C            IWORK - INTEGER
  135. C                    VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  136. C                    K ELEMENTS OF WHICH CONTAIN POINTERS
  137. C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS,
  138. C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
  139. C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
  140. C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND
  141. C                    K = LIMIT+1-LAST OTHERWISE
  142. C
  143. C            WORK  - DOUBLE PRECISION
  144. C                    VECTOR OF DIMENSION AT LEAST LENW
  145. C                    ON RETURN
  146. C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
  147. C                     END POINTS OF THE SUBINTERVALS IN THE
  148. C                     PARTITION OF (A,B),
  149. C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
  150. C                     THE RIGHT END POINTS,
  151. C                    WORK(LIMIT*2+1), ...,WORK(LIMIT*2+LAST) CONTAIN THE
  152. C                     INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
  153. C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3)
  154. C                     CONTAIN THE ERROR ESTIMATES.
  155. C***REFERENCES  (NONE)
  156. C***ROUTINES CALLED  DQAGIE,XERROR
  157. C***END PROLOGUE  DQAGI
  158. C
  159.       DOUBLE PRECISION ABSERR,BOUND,EPSABS,EPSREL,F,RESULT,WORK
  160.       INTEGER IER,INF,IWORK,LAST,LENW,LIMIT,LVL,L1,L2,L3,NEVAL
  161. C
  162.       DIMENSION IWORK(LIMIT),WORK(LENW)
  163. C
  164.       EXTERNAL F
  165. C
  166. C         CHECK VALIDITY OF LIMIT AND LENW.
  167. C
  168. C***FIRST EXECUTABLE STATEMENT  DQAGI
  169.       IER = 6
  170.       NEVAL = 0
  171.       LAST = 0
  172.       RESULT = 0.0D+00
  173.       ABSERR = 0.0D+00
  174.       IF(LIMIT.LT.1.OR.LENW.LT.LIMIT*4) GO TO 10
  175. C
  176. C         PREPARE CALL FOR DQAGIE.
  177. C
  178.       L1 = LIMIT+1
  179.       L2 = LIMIT+L1
  180.       L3 = LIMIT+L2
  181. C
  182.       CALL DQAGIE(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
  183.      *  NEVAL,IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST)
  184. C
  185. C         CALL ERROR HANDLER IF NECESSARY.
  186. C
  187.        LVL = 0
  188. 10    IF(IER.EQ.6) LVL = 1
  189.       IF(IER.GT.0) CALL XERROR(26HABNORMAL RETURN FROM DQAGI,26,IER,LVL)
  190.       RETURN
  191.       END
  192.